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The stability and causality of the Landau-Lifshitz theory and the Israel-Stewart type causal 
dissipative hydrodynamics are discussed. We show that the problem of acausality and instability 
are correlated in relativistic dissipative hydrodynamics and instability is induced by acausality. We 
further discuss the stability of the scaling solution. The scaling solution of the causal dissipative 
hydrodynamics can be unstable against inhomogeneous perturbations. 

I. INTRODUCTION 

f^ ' Presently the hydrodynamical approach is known as one of basic tools for the description of the collective aspects 

f^ of the relativistic heavy-ion collisions. However, instead of extensive analyses of experimental data based on this 
CN ' approach 1] , the study of the effect of dissipation is yet poorly explored [3, S lil • 

1— I [ Till now, two basically different approaches to describe relativistic dissipative hydrodynamics are commonly em- 

, ^ i ployed; one is the relativistic extension of Navier-Stokes theory introduced by Landau-Lifshitz (LL) [5] and Eckart 
Q, and the other is the causal dissipative (CD) hydrodynamics or second order theory, where the version due to 
^^ . Israel- Stewart Q is most popularly known. The relativistic Navier-Stokes theory (NS) is usually considered in the 
literature as a natural covariant generalization of the Navier-Stokes equation and in contrast to CD hydrodynamics, 
,— — 1 1 it is also refereed to as the first order theory. Authors who applies the first order theory argue that the second order 
,^ ■ theory is not yet completely established and the difference of two approaches would be not significant, since the effect 
Mh of viscosity is the measure of the deviation from the local thermodynamical equilibrium, and should be small where 
f^ ; the hydrodynamical approach is meaningful. However, a crucial point of the reason why the second order approach 
1^ ] should be explored is that the propagation speed in the first order theory is infinite and it does not satisfy the rela- 
te , tivistic causality. It is also known that the first order theory leads to dynamical instabilities. These aspects are not 
only the question of conceptions but also related to serious practical problems and deserve a more detailed analysis. 
I , One way to solve the problem of acausality in a first order theory is introduce, for example, a memory effect with 

^ ' a finite relaxation time [3[. In this way, we can derive a CD hydrodynamics from the LL theory. There are several 
C3 ] different approaches to derive the CD hydrodynamics 0,010, [9, 10, 11]. In this work, we consider the Israel-Stewart 
04 type CD hydrodynamics. 

^~~^ ' A crucial point is whether a relativistic dissipative hydrodynamics is stable or not. If a theory is unstable, it will 

■ , be very difficult to extract any physically meaningful results from it. The analysis of the stabilityof relativistic 

C^^ ■ dissipative hydrodynamics has been extensively studied by Hiscock and his collaborators [l2, [ij, [iJ, UM, \1m, ll3|- 

^5 , Their conclusions are summarized as follows. 1) The Eckart theory is unstable for the linear perturbation around 

the hydrostatic states whereas the LL theory and the IS theory (of the Landau frame) are stable [ij, [lj| . 2) The 

LL theory is shown to be unstable for the linear perturbation around hydrostatic states in a general frame where the 

fluid is Lorentz boosted [13| . 

Another important analysis was implemented by Kouno et al [l^- They discussed the linear perturbation around 
the scaling solution of the LL theory and found that the scaling solution of the LL theory can be unstable. A similar 
C^ ' discussion is repeated by [19| using another equation of state. 

The purpose of this paper is to complement these discussions. We will, in particular, focus on two subjects. One 
is the relation between causality and stability. In general, the concept of causality and stability are independent, but 
in relativistic systems, as we will see, they are closely related. To show this, we discuss the stability from Lorentz 
boosted frames and show that instability is induced because of acausality. 

The other is the stability around the scaling solution in the CD hydrodynamics. In the LL theory, the scaling 
solution can be unstable, although it is always stable when the Reynolds number is larger than one. The stability 
around the scaling solution has not yet been discussed in the CD hydrodynamics. The applicability of the scaling 
ansatz is not obvious in the CD hydrodynamics, because, as we will see later, the numerical calculation of the CD 
hydrodynamics shows a kind of non-periodic oscillations in the central rapidity region. 

In this paper, we restrict ourselves to the discussion on 1+1 dimensional motion of massless ideal gas (to be specific, 
an ideal three flavor massless QGP gas), for simplicity. This paper is organized as follows. In section[ni we discuss the 
stability and causality around the hydrostatic states. The result of this section has already been shown in [4|. I IS. 1 14| . 
To establish the relation between causality and stability, we discuss the stability of the hydrostatic state from Lorentz 
boosted frames in section Hill In section llVl the stability around the scahng solution is discussed. Section|V]is devoted 
to concluding remarks. 
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II. LINEAR PERTURBATION AROUND THE HYDROSTATIC STATE 

Before discussing the relation between causality and stability, we will discuss first the stability near the hydrostatic 
state following 0, [ll [13 . 

In the CD hydrodynamics of 1+1 dimensional systems, the equations are given by 

a.T^'' - 0, (1) 

Tfl-^n + n = -cd^u'^. (2) 

ar 
Here, 11 is the bulk viscosity and T'^'^ is the energy-momentum tensor defined by 

T^"" = {e + P + U)u^'u'' -{P + U)g'''', (3) 

where e, P and u^ are the energy density, pressure and four-velocity of the fluid, respectively. For the massless QGP, 
we have e = 3P, and e + P = Ts, where T and s are the temperature and entropy density, respectively. The fluid 
velocity is determined from the energy-momentum tensor following the definition of Landau-Lifshitz [5| . It should be 
noted that the CD hydrodynamics is reduced to the LL theory in the limit of the vanishing relaxation time r^j. 
For later convenience, we parametrize the velocity as follows, 

M^ = (cosh6',sinh6'). (4) 

We adopt the following parametrization of the bulk viscosity coefficient and the relaxation time [J] , 

C = as, (5) 

m - j^b, (6) 

where s is the entropy density. The parameters a and b characterize respectively the magnitudes of the viscosity 
and the relaxation time [23|. As we will see later, the parameter b should be smaller than 3/2 to be consistent with 
causality Q. 

We consider the plane wave perturbation around Eq, 6q and Hq, 

e{t,x) w so{t,x) +eei(t,x), (7) 

9{t,x) « 9o{t,x)+eei{t,x), (8) 

n{t,x) w Uo{t,x) + eUi{t,x), (9) 



El {t,x) 

6i{t,x) I =6^'^*-''=^ I ^1 I ■ (10) 




with 



ni(i,x) 

Here e is a small expansion parameter. The two transport coefficients are also expanded as 

C ~ Co + e<5C, (11) 

tr = ~ tro + e<5Tfl. (12) 

For the linear perturbation around the hydrostatic state, 

SOeqnarray* — COnSt, 6*0 = Hq = 0, (13) 

the evolution equation of the linear perturbation is given by 




where 



(14) 



iuj —ik{eo + Po) 
A^ I a{-ik) iLu{eo + Po) -ik I , (15) 



-ikQ 



.0 




where a = dP/de = 1/3. To have non-trivial solutions, the determinant of the matrix A should vanish so that the 
frequency uj has to satisfy the following dispersion relation, 

u;3-^u;^-(k^^+a]k^cj + i — k^ = 0. (16) 

tro \TReQ + Pq ) tro 

The behaviors of the frequency characterizes the stability and the propagation speed of the fluid. 

A. Landau-Lifshitz theory 

We consider the case of the LL theory by taking tr = 0. Then the solution of the dispersion relation p^ is 
analytically given by 



'^° ' ^fc^- ,, ^\,, k\ (17) 



2{eo + Po) V 4(e + Po) 

The real and imaginary parts of the frequency cu are shown in Fig. [T]for a — 0.1 and T = 200 MeV, respectively. 
One can easily see that the behavior of uj changes at the critical momentum kc — 1\fa(e + P)/Cq', below the critical 
momentum, there are two propagating modes, while they are changed to two non-propagating modes above it. 

We assume that the propagation speed of fluid is characterized by the group velocity for propagating modes. Then, 
the causality of the theory is determined by the behavior of the real parts of the frequencies. For the small k, the 
propagation speed is given by 

This is nothing but the usual sound velocity, and the LL theory seems to be consistent with causality. For the large 
k, however, the propagating modes are changed to the non-propagating modes which show k^ dependence. This 
momentum dependence is the same behavior as that of the non-propagating mode in diffusion processes, where the 
propagation speed is infinite. It is considered that the behavior of the non-propagating mode is the origin of acausality 
in the LL theory. 

On the other hand, the stability of the theory is characterized by the behaviors of the imaginary parts of the 
frequencies. One can easily see that the two modes always have positive imaginary parts, and hence the LL theory is 
stable under the linear perturbation around the hydrostatic states. 

Note that this is different behavior from the Eckart theory, where the theory is acausal and unstable even for the 
linear perturbation around the hydrostatic state [13| . 

B. Causal dissipative hydrodynamics 

Different from the case of the LL theory, we can still obtain the dispersion relation from E^. pB)) , but the analytic 
form of the solution becomes extremely complicated. However, in the large k limit, we have 



UJ = { V b ' " ' 2rn{l+ab) (I9) 

\ ■ ab ^ ^ 



k 'TRil+ab) 

while for the small k, 



UJ = < ^'^ 1-T^kyb "^ 2(h-T^fe2) _ (20) 

In this case, there are three modes; two of them are propagating modes and the remaining one is a non-propagating 
mode. From Eqs. p^ and (^0]) . we can see that, for the small k, the group velocity of the propagating modes reduces 
to that of the ideal one, -/a, like the LL theory. On the other hand, for the large k, the group velocity is given by 



V, = ^l/b + a. (21) 

That is, the group velocity is affected by the bulk viscosity. This gives the maximum propagation speed of this ffuid. 



- 0- 



0) 




100 



10 20 

k(1/fm) 



30 




10 20 

k(1/fm) 



30 



FIG. 1: The real (left panel) and imaginary (right panel) parts of the frequency in the Landau-Lifshitz theory at the rest frame. 
The two propagating modes (the solid and dotted lines) are changed to the non-propagating modes at the critical kc- Both of 
them have positive values. 
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FIG. 2: The real (left panel) and imaginary (right panel) parts of the frequency in the causal dissipative hydrodynamics at the 
rest frame for a = 0.1 and 6 = 6. There are three modes. One is non-propagating mode (the solid line), and the other two are 
propagating modes (the dotted lines). The two imaginary parts of the propagating modes are degenerated. 



In Fig. [2j we show the real and imaginary parts of the frequency w as functions of momentum fc for a = 0.1, 6 = 6, 
and the temperature T — 200 MeV. From the left panel, one can see that the group velocity of the two propagating 
modes dRe oj/dk are still slower than the speed of light. As for the non-propagating mode, we find, from the right 
panel, that the imaginary part becomes a constant for the large k, always remaining positive. This is true for the 
two propagating modes. The positivity of the imaginary part guarantees the stability of the hydrostatic state for 
plane-wave perturbations. That is, the CD hydrodynamics, with this parameter set, is causal and stable. 

However, as was mentioned, the propagation speed of the fluid (pij) is affected by the parameter b. For the ideal 
equation of state, a = 1/3, the propagation speed exceeds the speed of hght if we use the parameter b < 3/2. In 
Fig. [21 the real and imaginary parts of uj as functions of momentum k are plotted for a = 0.1,6 = 1, and the 
temperature T ~ 200Mey. In this acausal parameter set, Vc is larger than one, and hence, as one can see from the 
left panel, the propagation speed of the propagating modes dRe u/dk exceeds the speed of light. That is, even the 
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FIG. 3: The real (left panel) and imaginary (right panel) parts of the frequency in the causal dissipative hydrodynamics at the 
rest frame for a — 0.1 and 6=1. There are three modes. One is non-propagating mode (the solid line), and the other two are 
propagating modes (the dotted lines). The two imaginary parts of the propagating modes are degenerated. 

CD hydrodynamics can be acausal depending on parameter sets. However, all the modes have negative imaginary 
parts and the theory is still stable. 

From these results, one may consider that the problem of acausality is independent of that of instability. However, 
as we will see in the next section, both problems are correlated in relativistic systems. 

In this section, we discussed the propagation speed under the linear approximation. It should be noted that the 
propagation speed is changed when the non-linear effect is taken into account. See Appendix |X] for details. 

III. STABILITY IN GENERAL EQUILIBRIUM FRAME 

In the previous section, we discussed the stability of a small perturbation mode around the hydrostatic state. Then 
we found that even if the dynamics is not consistent with causality (the LL theory and the CD hydrodynamics with 
the acausal parameter set), the hydrostatic states are still stable. Thus stability is not related with the causality of the 
theory in these mode. However, these two should be related. Suppose an acausal propagation of a wave in a covariant 
theory is allowed. Then an initial pulse within the light-cone eventually would develop a singular behavior at the 
light-cone, since the light-cone cannot be crossed within a covariant theory. To clarify this point, we will investigate 
the behaviors of the perturbation in a general Lorentz boosted frame. 

Let us consider the linear perturbation around the hydrostatic state from the Lorentz boosted frame moving with 
the velocity V. Then the total velocity of the fluid is given by 



[/^ = 7v(cosh6'-|-l/sinh6', Fcosh6'-|-sinh( 
== cosh(V' + 6')(l,tanh(V' + 6')), 



(22) 



where tanh?/) = V. 

Substituting this into the energy-momentum tensor and repeating the same linear analysis around the hydrostatic 
state with ip being maintained constant, the evolution equation for the linear perturbation is given by 



£l 

A\ 9, 

Hi 



= 0, 



(23) 



where the components of the matrix are 

All ~ (cosh^ -0 + sinh^ -i/;a)(ia;) + cosh'0smh?/'(l + q;)(— i/c), (24) 

Ai2 — 2{eQ + Po)coshipsmh.^{iLu) + woicosh^ tp + sin]!^ ip){—ik), (25) 

Ai3 = sinh^ ipliu) + cos'h'ipsiiih'ip{—ik), (26) 

A21 — cosh ip sinh ■tjj {I + a) (iuj) + {sinh^ ip + cosh^ ip a) {—ik), (27) 

A22 — u'o(cosh^ V + sinh^ ■(/j)(ia;) + 2u'oCOsh?/;sinhV'(— ifc), (28) 

A23 — coshipsmh^{iu)) + cosh^ i/j{—ik), (29) 

A31 = 0, (30) 

A32 — Co(sinh-0(Jti;) + cosh'!/;(— iA:)), (31) 

^33 = tro coshipiiuj) + tro sinhip{-ik) + 1. (32) 

Then the dispersion relation is obtained by solving the following equation, 

iAuj^ + iBkuj^ + Clu^ + iDk^uj + Ekuj + iFk^ + Gk^ = 0, (33) 

where 

A = cosh6l(-l- (l-w2)gjj^jj2^.|^ ^.g^^l 

B = sinhV'(l-(l-Wc) + 3(l-Wc)cosh2V), (35) 

C = ^(a + (l-a)cosh2V), (36) 

tro 

D = coshV'(-3(l-w2)cosh2V + 2(l-Wc) + l), (37) 

2 

E == (l-a)cosh-0sinhV', (38) 

TRO 

i^ = sinhVX-l + (l-Wc)cosh^V), (39) 

G = (1- (l-a)cosh2V)- (40) 

Tro 

One can easily see that the behavior of the frequency uj depends on the choice of Vc, which is the propagation speed 
of the fluid defined by Eq. (^1]) . 

A. Landau-Lifshitz theory 

As we discussed, the Landau-Lifshitz theory has two modes in the local rest frame. In the Lorentz boosted frame, 

however, we have three modes, because of the following reason. In this case, the coefficients of Ec^. ([55]) are given by; 

A = CoCOshV'sinh^i/), (41) 

B = CosinhV'(l-3cosh2V'), (42) 

C = -ieo + Po){a + {l-a)cosh^i;), (43) 

D = Cocosh-0(-2 + 3cosh^V), (44) 

E ^ 2(eo + ^o)(l -a)cosh7/'sinhV', (45) 

F = — Cosinh-^cosh^V: (46) 

G ^ (£o + Po)(l-(l-a)cosh2v). (47) 

The coefficient A disappears only in the rest frame {ip ^ 0). That is, there exists a gap in the calculations of the rest 
frame and the moving frame. 

In Fig. m the real and imaginary parts of the frequency u! are shown for a = 0.1 and T = 200 MeV at V — 0.1. 
One can see that all the three modes are propagating modes. The real parts of the two propagating modes denoted 

by the dotted line are degenerated, but the imaginary part of one of them is negative. Thus the LL theory is unstable 
in the Lorentz boosted frame. 
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FIG. 4: The real (left panel) and imaginary (right panel) parts of the frequency in the Landau-Lifshitz theory in the Lorentz 
boosted frame with the velocity v — 0.1. There are three propagating modes. One of the modes is denoted by the solid line. 
The remaining two modes are degenerated and plotted by the dotted line. One of the imaginary parts (one of the dotted lines) 
is negative. 
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FIG. 5: The real ( left panel) and imaginary ( right panel) parts of the frequency in the causal dissipative hydrodynamics in 
the Lorentz boosted frame with the velocity v = 0.9 for a — 0.1 and b — 1. There are three modes denoted by the solid, dotted 
and dashed lines, respectively. All modes have positive imaginary parts. 



B. Causal dissipative hydrodynamics 



First, we consider the causal parameter set, a — 0.1 and 6 = 6 where the propagation speed (|2ip is slower than the 
speed of light. In Fig. [5l the real and imaginary parts of the frequency u> are shown for T = 200 MeV at F = 0.9. 
From the behaviors of the real parts, one can see that the three group velocities become close to the speed of light, but 
never exceed. On the other hand, all three propagating modes have positive imaginary parts. As far as we checked, 
this theory does not becomes acausal and does not have a negative imaginary part for any V. Thus in this causal 
parameter set, the CD hydrodynamics is consistent with causality and stable even in the Lorentz boosted frame. 

However, this is not true for the acausal parameter set, for example, a = 0.1 and b = 1. In Fig. [51 we show the real 
and imaginary parts of the frequency lo. There are three propagating modes denoted by the solid, dashed and dotted 
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FIG. 6: The real (left panel) and imaginary (right panel) parts of the frequency in the causal dissipative hydrodynamics in the 
Lorentz boosted frame with the velocity v = 0.9 for a = 1 and b = 1. There are three modes denoted by the solid, dotted and 
dashed lines, respectively. One mode (the dotted line) has a negative imaginary part. 

line. It is clear that the group velocity is faster than the speed of light. Interestingly enough, one of the propagating 
modes denoted by the dotted line has a negative imaginary part. Thus the CD hydrodynamics with acausal parameter 
set is unstable in the Lorentz boosted frame. 

This result means that causality and stability are correlated and instability is induced by acausality in the relativistic 
dissipative hydrodynamics. A consistent theory should not change its stability depending on the choice of the frames. 
Thus, the LL theory is not consistent theories as candidates for the relativistic dissipative hydrodynamics. This is so 
also for the CD hydrodynamics with acausal parameter sets. 

We should stress, as a matter of fact, that we cannot implement stable numerical calculations of the CD hydrody- 
namics when we use acausal parameter sets. 

The stability of the LL theory from a general frame is discussed also in [13!| in a different context. See Appendix [Bl 
for details. 

IV. STABILITY AROUND SCALING SOLUTION 

In this section, we discuss the stability around the scaling solution. The scaling variables r and y are defined by 



= ^f^ - Z2 



By using these variables, the equations for the conservation of energy-momentum are reexpressed as 





1 , 


\t + z^ 


1/ = 


= -hl 






2 


[t- z\ 



{rdr + tanh(6' -y± Bs)dy){(t> ± B) 



cosh(6' ~ y)±Cs sinh(6' — y) 



Rr^ 



c,ve ±i De + -vn ) 



= 0, 



(48) 



(49) 



where tanh6's = \/a- Here we define the Reynolds number R ^ = — 11/ (Ts), which reproduces the definition of [16 
in the vanishing relaxation time limit. The new variable (j) satisfies the following relation, 



^/adlns = —=d\TiT. 



(50) 



The equation for the viscosity is 



Tndr^ + n = -C sinh(6l - y)dr + cosh(e' - y)-dy 9 



(51) 



We assume that the velocity of the fluid is given by 



tanh^ = — = tanhw. 
t " 



(52) 



This scahng ansatz is considered to be vahd near the central rapidity region. Thus, the equations of the scaling 
solution are given by setting 9 = y, 



Tdr(j)0 + (1 - Ro')c" - 0, 
T 



TBodrUa + Uq 



(53) 
(54) 



In the LL limit {tji -^ 0) , the equations are reduced to those obtained in [1^ . 

To see the stability of the scaling solution, we consider the linear perturbation as follows. 



9 = y + S9, (55) 

<j> = <l)o + S(l>. (56) 

Substituting them into Eqs. (|49p and (|5ip . the evolution equations of the linear perturbations are given by 

'dc. 



TdM + (1 - Ro^)c°dyS9 + {1~ R^^) ( - 



SR-'c';, = 0, 



Rn 



1 Rn 



r(l - R^')drS9 + c^.dyS^ + (1 - (cO)2)^0 - (1 - ic",r)S9Ra' - -^dySU + — ^^(tHq + Co 



no Tn 



-«'- + -H-^)(n.4«-^(^)^..-C.^^ 



where 



SR-^ = - 



1 



-sn 



d(t> Jo 
Ho fdln{Ts) 



(Ts)o iTs)o V d^ 



with the scaling solutions (j>o and Ho. 
By using the Fourier transform for y, 



A(r,y) = / dke-'''yA{T,k), 



Eqs. (|57|) . (I58|l and ([59)) are summarized as 



where 



TdrX = AX, 



X = 



5(f) 



(Jinn 



and 



/-(l-V)(^)-^o '(! + ") -ii-Ro')V^i-^k) V^R^' \ 



A 



(-»fc) 



V 



where f = t/tj^q and 






91n(rs) 



-(i-a)- ^^::,^^ i?,-^^^ 



ik 1 



d(j) 
SlnC 



(57) 

0, (58) 
(59) 

(60) 



(61) 



(62) 



(63) 



(64) 



/a 



(65) 
(66) 
(67) 
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FIG. 7: The phase diagram of the stabihty in the LL 
theory as a function of k and Ro- There are two stable 
regions. 



FIG. 8; The phase diagram of the stabihty in the CD 
hydrodynamics as a function of f and Ro- We set k — 
and 6 = 6. There are two stable regions. Even though 
Ro > 1, the scaling solution of the CD hydrodynamics 
can be unstable. 



Note that we consider the ideal gas equation of state where a = 1/3 and {d^/a/dcj))^ = 0. 

One can see that the scaling solution is unstable when f > 1/fe, for fc = and Rq = 1 + e, and r < 1/&, for 
fc = and Rq = 1 — e, because the equation for 59 becomes decoupled and A22 becomes positive. Here e is a small 
arbitrary constant. 

In general, the analysis of the stability is non-trivial unless the matrix A can be diagonalized. We discuss the stability 
in the Lyapunov direct method used in ^IS] . In this approach, the stability of the theory is analyzed by introducing 
the Lyapunov function, which characterizes the deviation from the scaling solution. The Lyapunov function should 
be 1) positive definite and 2) monotonically decreasing function with respect to the measure of the distance of the 
perturbed solution from the non-perturbed one. If we can find the Lyapunov function at given fc, 6, f and Rq, the 
scaling solution is stable for the parameter set. If the assumed Lyapunov function is found to be a monotonically 
increasing function then the scaling solution is unstable for the respective parameter set. See Appendix [C] for details. 

However, we should be noted that the stable region predicted in the Lyapunov direct method is normally under- 
estimated, unless we exhaust every possible Lyapunov functions. As an example, let us consider the limit of the LL 
theory which is realized when we set tr = dlnTu/dcj) = 0, in Eqs. ((57|) . ([58| and ([59|l . In Fig. [71 we show the phase 
diagram for the stability, on the k — Rq^ plane, assuming the Lyapunov function V = \6<j)\'^ + \S9\^. There are two 
stable regions; one is in i?o > 1 and the other in the small k and small Rq- We found that the scaling solution is 
stable also on the line of _Ro = 1 by solving the same V , but it is not shown in the phase diagram for simplicity. To 
confirm more precisely the stable regions, we have to study various possible Lyapunov functions. Then the real stable 
regions, in general, can distribute in broader regions on the phase diagram. As a matter of fact, the phase diagram 
obtained by Kouno et al. shows that the LL theory is stable for any k in all region of i?o ^ 1 [3 ■ On the other hand, 
we could not find unstable regions for this Lyapunov function. The stability of the domain between the two stable 
regions is not confirmed. However, from the analysis of Kouno et al., a part of the unconfirmed region should be an 
unstable region [ISJ . 

To analize the CD hydrodynamics, we first consider the following Lyapunov function V, 



V = X'fX = {6<pf + [Sef + {SlnUf 



Note that the time evolution of V is given by. 



TdrV = X''{A'' +A)X. 



(68) 



(69) 



For Rq = 1, we can eliminate the variable S(j) from the equations. In this case, we will consider the simpler Lyapunov 
function, 



V = Y^Y= {S(t)f + {6\nU)' 



(70) 
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FIG. 9: The phase diagram of the stabihty in the CD 
hydrodynamics which is calculated with the Lyapunov 
function V'" . We set k — and & = 6. Compared to Fig. 
([8|). the stable region is enlarged. 



FIG. 10: The phase diagram of the stability in the CD 
hydrodynamics which is calculated with the Lyapunov 
function V'" . We set fc = 1 and 6 = 6. The stable 
region for homogeneous perturbation (fc = 0) at low f 
around Ro = 1 is changed to the unconfirmed region for 
the inhomogeneous perturbation (fc = 1). 



Here the evolution equation of Y is given by 



where 



TdrY = BY, 



Y 



(5hin 



and 



B = 



-{1 + a) 

b) ^ v^fc 



bf-1 



(71) 



(72) 



(73) 



fcf-i 



In this case, we should discuss the eigen values of B^ + B. 
There are three parameters, k, f and Ro, fixing 6 = 6. 



First, we discuss the homogeneous perturbation setting 



k = and calculate the phase diagram in the f — Rq plane, as is shown in Fig. [H There are two stable regions; one 
is very small and located in i?o < I7 and the other is larger and located in _Ro > 1- This result is consistent with the 
instability of the scaling solution in f > 1/6 for Rq^ = 1 + e and f < 1/6 for R^^ = 1 — e. It should, however, be 
noted that one can see that the scaling solution is always stable on Rf^^ = 0, which can be seen from the behavior 
of Eq. (|5^ itself. On the other hand, we could not find unstable regions. As a matter of fact, we will discuss the 
stability by using various Lyapunov functions in the following, but still cannot find any unstable regions. 

As we have pointed out, the Lyapunov direct method normally underestimates the stable region. To fix the 
stable region, we have to discuss as many Lyapunov functions as possible. In this work, we analyzed the phase 
diagram with more three different functions: V = |(51ns|2 + |(56'|2 + |,51nn|2, V" = \S\ns\^ + \Se\^ + |i?fl ^^InHp and 
V" — \S(f>\'^ + \S9\'^ + IRq^SIuH]^. We found that the qualitative structure of the phase diagram is independent of the 
choice of these functions. We show the result of V'" in Fig. [9l for which we obtained the maximum stable region (See 
Appendix |D] for detailed form of the equation) . This phase diagram shows that most part of the phase diagram in 
i?o < 1 still belongs to the unconfirmed region. As for the region oi Rq > 1, we found that the stable region strongly 
depends on the choice of the Lyapunov function, and most of the physically accessible region is confirmed to be stable. 
As a matter of fact, the dashed lines in Fig. [9] shows the trajectories of the scaling solutions for a = 0.1 (left) and 
a = 1 (right), and one can see that most of the trajectories belongs to the stable region in the phase diagram. 

In this sense, we conclude that the scaling solution is stable for homogeneous perturbation. To see the stability for 
the inhomogeneous perturbation, we need to discuss the phase diagram for finite k. In Fig. 1101 we show the phase 
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FIG. 11: The numerical simulation of the 1+1 dimen- FIG. 12: The numerical simulation of the 1+1 dimen- 
sional CD hydrodynamics with a = 1 for t = 0.52 0.72 sional CD hydrodynamics with a = 1 for t — 0.52 0.72 
and 0.92 fm. The non-periodic oscillations appear in the and 0.92 fm. The non-periodic oscillations appear in the 
center of the fluid. center of the fluid. 

diagram for k = 1, which is calculated by using the function V'" . One can see that the stable region in the low f, 
near i?o = 1, is changed into a unconfirmed region. This behavior is commonly seen in the results obtained using 
other three Lyapunov functions. As k increases, this propensity becomes prominent and the stable region shrinks 
increasingly. The trajectories of the scaling solution are plotted in the same figure, again. One can see that the 
trajectory of the scaling solution with a = 1 passes the unconfirmed region, although the trajectory with a = 0.1 
still stays inside the stable region. It should be noted that, for A: = 0, at least a part of the unconfirmed region near 
Rq = 1 should be unstable. This suggests that the scaling solution can become more unstable for the inhomogeneous 
perturbation as the bulk viscosity a increases. 

The situation discussed here seems to be realized in the numerical simulations of the 1+1 dimensional CD hydrody- 
namics. When we increases the bulk viscosity a, we found that the numerical calculation becomes unstable and a kind 
of non-periodic oscillation appears in the center of the fiuid. In Figs. [TTjandfT^ we show the temperature and the 
bulk viscosity of the viscous fluid with a = 1 for i = 0.52 0.72 and 0.92 fm, respectively. We used the Landau initial 
condition with the initial temperature T = 590 MeV and the initial size 0.7 fm. To remove numerical oscillations 
which will disappear in the continuous limit, we used the additional viscosity which is introduced in [4]. It seems that 
the oscillation appears when Rq exceeds a critical value by decreasing the temperature and the bulk viscosity. The 
amplitude of the oscillation grows up with time and finally the numerical calculation collapses. 

So far, we discussed fixing the parameter b — 6 which is consistent with causality. Even when we use the acausal 
parameter set, 6=1, quantitative behavior of the phase diagram is not changed, but the stable region becomes much 
smaller. 



CONCLUDING REMARKS 



Although the physical importance is recognized for the application in QGP physics, the relativistic viscous hydro- 
dynamics is not well established yet. Some authors use the first order theory to estimate the effect of viscosity in 
collective observables such as V2 hoping that the deviation from the ideal hydrodynamics is small so that the theory is 
of the first order or the second order might be irrelevant. However, the difference between them might be fatal when 
any instabilities or singularities emerge, such as shock wave propagations. Therefore, it is fundamental to understand 
the stability of these theories. In this paper, we discussed the causality and stability of the two cases of relativistic 
dissipative hydrodynamics, the LL theory and CD hydrodynamics. 
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hydrostatic state 

moving frame 

scaling solution 



LL theory (acausal) 



stable 

unstable 

stable/unstable 



CD hydrodynamics (acausal) 



stable 

unstable 

stable/unstable 



CD hydrodynamics (causal) 



stable 

stable 

stable/unstable 



The LL theory is known to be acausal whereas the CD hydrodynamics can be causal depending on the values of 
parameters of the theory. The stability of the theories are summarized in the above Table. Around the hydrostatic 
state, the LL theory and the CD hydrodynamics are stable. However, when we move to a Lorentz boosted frame, 
the acausal theories like the LL theory and the CD hydrodynamics with acausal parameter set become unstable. The 
second line shows that causality and stability are intimately correlated in relativistic dissipative hydrodynamics. The 
stability of a theory should not depend on the choice of frames. In this sense, the LL theory and the CD hydrodynamics 
with acausal parameter sets are inconsistent. 

The stability of the scaling solution was analyzed by using the Lyapunov direct method. In the LL theory, it is 
known that the scaling solution is stable against homogeneous and inhomogeneous perturbations when we use initial 
conditions which satisfies i?o > 1 [l8|. In the CD hydrodynamics, we found that the scaling solution cannot be stable 
even for i?o > 1. For the homogeneous perturbation (fc = 0), we confirmed that most parts of the trajectories of the 
scaling solutions pass the stable region in the phase diagram, which was plotted in terms of f and Rq . Thus, the 
scaling solution will be stable for the homogeneous perturbation. However, as the k increases, the confirmed stable 
region in the phase diagram shrinks and the trajectories of the scaling solutions start to penetrate the unconfirmed 
region. When the unconfirmed regions are real unstable regions, it means that the scaling solution is unstable for 
inhomogeneous perturbations. This instability is distinguished for larger bulk viscosity because the trajectory with 
larger viscosity is easier to penetrate the unstable region. 

Above conclusion may be supported by the numerical calculations. As a matter of fact, we found that the numerical 
calculations of the 1+1 dimensional CD hydrodynamics becomes unstable as the bulk viscosity coefficient increases 
and a kind of non-periodic oscillations appears in the central rapidity region. To see the quantitative signature of the 
oscillations, we have to investigate various cases with different parameters and initial conditions. This oscillation may 
be interpreted as turbulence because the instability of the scaling solution indicates chaos which acts as the trigger 
of turbulence. However, in this work, we could not confirm the unstable region on the phase diagram. To see the 
appearance of turbulence, we need more systematic study of the instability around the scaling solution beyond the 
Lyapunov direct method. This is a challenge for the future. 

When we discuss the shear viscosity, we have to find the parameters which is consisitent with causality as was 
discussed in this paper. This is also a future task. 
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APPENDIX A: NON-LINEAR EFFECT FOR CAUSALITY 



The propagation speed of the fluid has been discussed based on the linear analysis, and hence the effect of nonlin- 
earity is ignored. In this section, following the discussion of [la, U^, we derive the effect of the nonlinearity in the 
propagation speed. 

For the simple 1+1 dimensional system, the hydrodynamic equations are summarized as follows; 



where y^ — (£,^,H) and 

(AY 



{A'^ydtY'' + [A^fd^Y" + B^ = 0, 



cosh^ + (cosh^ 9 -l)a 2w cosh sinh 9 sinh^ 9 

cosh^sinh6'(l + a) z/;(cosh + sinh 9) cosh^sinh^ 
C sinh 9 tr cosh 9 



(Al) 



(A2) 



[AY 



' cosh sinh 0(1 + a) w(cosh + sinh 9) cosh sinh 6* 



sinh" ( 



(sinh 9 + l)a 2w cosh 6* sinh 



V 







C cosh i 



cosh^ 9 
Tji sinh 9 



(A3) 
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5'"== (0,0,0, n). (A4) 

Here, w is the effective enthalpy density, w = e + P + 11. 
The characteristic speed v is given by 

det{viAy - (A)^) = 0. (A5) 

The speed is easily estimated in the local rest frame, 9 = 0. Then we have the following three solutions. 






One can easily see that if 11 is small and we can replace w with e + P, this result is same as Eq. ([T9| . 

We consider the case of the effective enthalpy density is positive. Then, to satisfy causality, the transport coefficients 
should satisfy the following condition, 

— <{l~a){e + P + U). (A7) 

TB. 

This is, again, the generalization of the restriction for the transport coefScients discussed below Eq. (|5|). 

APPENDIX B: INSTABILITY IN GENERAL EQUILIBRIUM FRAME (HISCOCK-LINDBLAM) 

In section , the stability from a Lorentz boosted frame was discussed. A similar problem was discussed by Hiscock 
and Lindblam [ij]. In this appendix, we apply their discussion to the CD hydrodynamics. 

They consider the transformation of the coordinate by using the following replacement of the variables, 

u) = -f{oJ + vk), (Bl) 

k = jik + vCj), (B2) 

where uj and k are variables in the rest frame, and lu and k are in the new frame, which moves with the velocity v. 
Substituting into the result obtained in the rest frame (J16p . we have 



C 1 , ^\ .3n„ , ~^2r~. , ,.l\ , „• Q^ ..2/K , ~\2 



rii^ + vky 7^(J) + vky - + a r {k + vQYiQ + vk) + i—Y{k + v(^Y = 0. (B3) 

TR \TRe + P J TR 

We can easily solve the equation for fc = 0, 

i 1 — OiiP' 

The imaginary part is positive and hence the theory is still stable. 
On the other hand, in the LL theory, the solutions are given by 

, , {e + P){l-av^) 

^ = 0,0,-1 —^ . (B5) 



Thus, the LL theory is unstable again. 



APPENDIX C: LYAPUNOV DIRECT METHOD 

Here, we summarize the stability analysis based on the Lyapunov function. As an example, let us consider the 
damped harmonic oscillator, 

|. = -- (CI) 

d 

— V = —JV — LJ^X. (C2) 
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The solution of the equation converges to a; = t; = 0. 

To discuss the stabihty around the equihbrium solution xq = Wq = 0, we introduce a function which characterizes 
the deviation from the equilibrium. For example, we choose 



V ^ [v-vof +a^{x~XQf. 



(C3) 



This is positive definite and if this function monotonically decreases with time, the system is stable and the function 
V is called the Lyapunov function. The time evolution of the function V is given by 



— V = {a{x~xo),v-vo)]\/n ^ "' 



dt 



V- uo 



where 



M = 



-27 



The eigen values of the matrix M are given by 



A± = -(-7± V7^ + (a-wVa)2). 



(C4) 



(C5) 



(C6) 



One can see that when a = w, T^ is a monotonically decreasing function in time. Thus V is the Lyapunov function 
and the equilibrium state is stable. 

However, if we use a ^ uj, the function generally have a positive and negative eigen values. Thus we cannot 
determine the stability of the equilibrium solution. In this sense, the Lyapunov direct method usually underestimates 
the stability of the system. 

Similarly, when we find that the minimum eigen value is positive and hence T^ is a monotonically increasing function, 
the equilibrium solution is unstable. 



APPENDIX D: ANOTHER CASE OF THE FUNCTION V" 



Instead of Eq. (|63p . we introduce the following vector. 



X 



R^^6\nIV 



Then, the Lyapunov function is given by 

Then the evolution equation of the Lyapunov function is 

TdrV" =X'^{A'^ +A)X, 

where 

/ -(1 - R^') (%f )^ - R^' (1 + a) -(1 - R-^)^{-^k) 



A = 



I— i-ik) 



—ik 



Similarly, as for Rq = I, the matrix B is given by 

-(1 + a) 
1 



B = 



^ir-l) + -h + ^ -- ' 



/a 

(-»fc) 



-T - [1 - R-')il + a) J 



(Dl) 

(D2) 
(D3) 

(D4) 



(D5) 



fcf-i 
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